#ifndef DIST_POLYA_HEADER
#define DIST_POLYA_HEADER

#include<string>
#include<fstream>
#include<sstream>
#include<iostream>
#include<algorithm>
#include<iterator>
#include<vector>
#include<map>
#include<cstdlib>
#include<ctime>
#include<cmath>
#include <gsl/gsl_sf_gamma.h>
#include "def.hpp"
#include "utils.hpp"

using namespace std;

double polya_likelihood(int**& data, int*& n_i, double*& param, int& K, int& N){
  double paramSum = get_sum(param, K);
  double likelihood = N * gsl_sf_lngamma(paramSum);
  for(int i=0; i < N; i++){
    likelihood = likelihood - gsl_sf_lngamma(n_i[i] + paramSum);
    for(int k=0; k < K; k++){
      likelihood = likelihood + gsl_sf_lngamma(data[i][k] + param[k]);
      likelihood = likelihood - gsl_sf_lngamma(param[k]);
    }
  }
  return likelihood;
}

double polya_likelihood(int**& data, int*& n_i, double*& param, int& K, int& N, int*& id_list){
  double paramSum = get_sum(param, K);
  double likelihood = N * gsl_sf_lngamma(paramSum);
  for(int i=0; i < N; i++){
    likelihood = likelihood - gsl_sf_lngamma(n_i[id_list[i]] + paramSum);
    for(int k=0; k < K; k++){
      likelihood = likelihood + gsl_sf_lngamma(data[id_list[i]][k] + param[k]);
      likelihood = likelihood - gsl_sf_lngamma(param[k]);
    }
  }
  return likelihood;
}


#endif
